Anomalous structural evolution of soft particles: Equibrium liquid state theory 
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We use the hyper-netted chain approximation of liquid state theory to analyze the evolution with 
density of the pair correlation function in a model of soft spheres with harmonic repulsion. As 
observed in recent experiments on jammed soft particles, theory predicts an 'anomalous' (nonmono- 
tonic) evolution of the intensity of the first peak when density is increased at constant temperature. 
This structural anomaly is a direct consequence of particle softness, and can be explained from 
purely equilibrium considerations, emphasizing the generality of the phenomenon. This anomaly is 
also predicted to have a non-trivial, 'S-shaped', evolution with temperature, as a result of a com- 
petition between three distinct effects, which we describe in detail. Computer simulations support 
our predictions. 

PACS numbers: 05.20.Jj, 64.70.qd, 83.80.Fg 



In very recent work 0, 0, H, 3, the structure of soft 
repulsive particles was analyzed experimentally and nu- 
merically. Particular attention was paid to a nonmono- 
tonic behaviour of the first peak of the pair correlation 
with increasing density, as shown in Fig. [T] In the limit 
where particles arc infinitely hard, the height of the first 
peak diverges when approaching a maximal packing frac- 
tion where the pressure diverges [B|, 0] , and configurations 
with larger density cannot be explored. Thus, the non- 
monotonic height of the pair correlation function g(r) was 
understood as a smooth, finite-temperature, signature of 
this jamming transition [lj. While simulations study soft 
repulsive potentials at finite temperatures [l|, Q , experi- 
ments introduce instead particle softness, which similarly 
allows for finite overlaps between particles. For instance, 
in the very elegant experiment of Ref. 0], the pair cor- 
relation function of a slowly swelling assembly of (com- 
pressible) tapioca pearls is followed, see Fig. [TJ 

These studies were performed in noncquilibrium condi- 
tions following a preparation protocol which varies from 
one experiment to another. Here, 'noncquilibrium' means 
that particles rearrangements are rare during the course 
of an experiment, and do not allow the system to sample 
equilibrium trajectories at any given state point, even 
at finite temperature. The jamming transition cannot 
be probed at thermal equilibrium, because the system 
first hits a glass transition and cannot maintain ergodic 
behaviour at such large densities @, @, @]- Thus, the 
structural evolution along particular noncquilibrium tra- 
jectories does not result from ensemble averages in the 
usual sense. When comparing two different state points, 
it is understood that a single configuration is transported 
from one state point to another along a specific path. 
This experimental procedure suggests that sample prepa- 
ration protocols might well play an important role. 

There are several questions left unanswered by these 



studies. First, the origin of the maximum is systemati- 
cally attributed to an underlying T = jamming transi- 
tion, and it is thus unclear how the phenomenon should 
be described theoretically. Second, the effects of sample 
preparation were either not discussed, or claimed to be 
negligible [l|. This seems to contradict recent reports 
that the location of the jam ming transition is sensitively 
protocol dependent [n], II, fl2l[l3| a conclusion which 



should extend to finite temperatures as well. Third, the 
location of the maximum was predicted to obey a simple 
scaling behaviour near zero temperature [l[ , as confirmed 
by simulations spanning a large temperature window 0] • 
However, this prediction must break down when temper- 
ature is large enough to allow relaxation in the system, 
but this has not been explored. 




FIG. 1: Nonmonotonic evolution with increasing density of 
the height of the first peak in the pair correlation of a slowly 
swelling assembly of (soft) tapioca pearls. Volume fractions tp 
are normalized with respect to that of the maximum location 
at (fi c . The inset shows the height of the first peak near r « 1.1 
as a function of <p/ip c - These data are adapted from Ref. [3], 
and were kindly provided by X. Cheng. 
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In this paper, our aim is to gain theoretical insight 
into the behaviour of the pair correlation function of 
dense assemblies of soft repulsive particles and to ad- 
dress the three above questions. We suggest that the 
simplest starting point is equilibrium statistical mechan- 
ics, and we use the tools of liquid state theory to attack 
the problem. Although we know beforehand that our re- 
sults should eventually become inapplicable to the low 
temperature regime where simulations and experiments 
fail to probe equilibrium configurations, our simple ap- 
proach provides useful insights, which we report in this 
short note. 

We follow Refs. [H, El and study an assembly of N soft 
repulsive particles interacting through harmonic repul- 
sion: 



V(r l} < 1) = e(l - r, y -) 2 



(1) 



where = \vi — rj 



denotes the distance between par- 
ticles i and j. Particles separated by > 1 do not 
interact, V(tyj > 1) = 0. In the following, we report 
temperatures in units of e, and the particle diameter sets 
the unit lengthscale. We study the pair correlation func- 
tion defined by 
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where the brackets indicate an ensemble average at ther- 
mal equilibrium, and p = N/V is the number density. 
We use the hyper-netted chain (HNC) closure relation 
for the pair correlation function [14| . 



?(r) = exp[— pV(r) + g(r) - 1 - c(r)], 



(3) 



for the pair potential V(r) in Eq. J!}. We have de- 
fined (3 = 1/T, where T is the temperature. In Eq. ([3]), 
c(r) is the direct correlation function defined through the 
Ornstein-Zernike equation: 

g(r) - 1 - c(r) +pf dr'c(r - r')O(r') - 1). (4) 

Although HNC is known to become quantitatively in- 
accurate in the hard sphere limit at the large densities 
investigated here [3, it is a convenient tool to study 
the interplay between temperature and dynamics in soft 
repulsive spheres. For the relatively large temperatures 
studied below, we shall find that HNC in fact compares 
quite well with simulations, at least qualitatively. 

We solve the HNC closure relation in Eq. © numer- 
ically. To this end, we discretize g(r) between r = 
and r = 32, using 2 15 points and use a standard Pi- 
card iterative procedure, as described in more details in 
Ref. [HI]. As a result, we obtain g(r) over a broad range 
of temperatures, T € [0.00076,0.5] and volume fractions, 
p = np/6e [0.5,1.2]. 
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FIG. 2: Nonmonotonic evolution of the height of the first peak 
in the pair correlation g(r), obtained using the hyper-netted 
chain approximation. The temperature T is fixed, and the 
volume fraction increases from right to left with ip = 0.600, 
0.663, 0.777, 0.900, and 1.00. 



Our first result is shown in Fig. [21 where the pair cor- 
relation of harmonic spheres at constant temperature, 
T = 1.5 • 10 -3 , is shown for volume fractions increasing 
from ip = 0.60 up to ip — 1.0. While the position of the 
first peak shifts to smaller distances as a result of the 
compression, the height of the first peak displays a non- 
monotonic evolution. For this temperature, it increases 
up to tp w 0.777, and then decreases for larger volume 
fractions. This evolution is very similar to the one re- 
ported in Fig. [TJ In Fig. [2j it can be seen that the sec- 
ond peak of g(r) also displays a nonmonotonic behaviour, 
while we find the subsequent peaks do not. Remarkably, 
the same observations were reported in experiments Q. 

In the context of jamming and glassy physics, this re- 
sult is unexpected, since in simple glass-forming liquids, 
modelled for instance as Lennard-Jones fluids, the first 
peak of g{r) typically increases with density at large den- 
sity, contrary to the system of harmonic spheres, where 
it instead decreases, see Fig. O 

This result shows that a nonmonotonic structural evo- 
lution of g(r) exists even at thermal equilibrium and not 
only in the nonequilibrium iso-configurational conditions 
investigated in Refs. [3, [2, 0, 9 . Therefore, we conclude 
that equilibrium concepts can be used to understand its 
origin and we suggest that the mechanism behind the 
behaviour in Fig. [2] is quite generic and not specifically 
due to an underlying nonequilibrium jamming transition 
at T = 0. Indeed, we find in the literature several ob- 
servations of similar anomalies observed in equilibrium 
conditions in systems of soft repulsive particles, for in- 
stance for Hertzian spheres [III or in the Gaussian core 
model 
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Il8j . To our knowledge, the close correspon- 
dence between noncquilibrium observations near the jam- 
ming transition, and these earlier reports at thermal equi- 
librium was not noted before. 

Physically, the anomaly is a direct consequence of par- 
ticle softness 0, 17, 18, [l9|. The system behaves as 



3 



io- 



10 



-2 



E-H 



10 



-3 _ 



10 



-4 



~i 



0.68 



0.72 



0.76 



0.8 



FIG. 3: Location of the maximum of the first peak of g(r) in 
a temperature (T) - volume fraction (ip) phase diagram. This 
location is itself a nonmonotonic function of temperature. In 
particular in the low T regime, it shifts to large ip when T 
decreases, as the system is able to find more efficient ways to 
pack non-overlapping particles. 



an effective hard sphere fluid at sufficiently low volume 
fraction where it can easily minimize energy by avoid- 
ing particle overlaps. In this regime, the height of the 
peak increases with ip, just as it does in a hard sphere 
fluid. When the volume fraction increases, however, it 
becomes harder, and thus entropically less favourable, 
to avoid particle overlaps as the system must find more 
efficient ways to pack the particles. Thus, the system 
may gain free energy by increasing the disorder at the 
expense of allowing particle overlaps, which are energet- 
ically not very costly for soft repulsions such as the one 
in Eq. ((T|). In this regime, disorder increases when in- 
creasing ip, and this broadens the first peak of g(r) and 
decreases its height. We note that similar competitions 
are invoked to explain an even larger set of anomalous 
thermodynamic behaviours, as observed for instance in 
widely studied network-forming liquids such as silica [20I ] 
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and water 

Since we attribute the anomaly to a competition be- 
tween energy and entropy, the outcome of this compe- 
tition can be expected to have an interesting evolution 
with temperature. Thus, we now turn to the influence of 
temperature on the nonmonotonic behaviour described 
in Fig. [2 In Fig. [3J we follow the locus of the maxi- 
mum of g(r) in a (T, ip) phase diagram, as predicted from 
the HNC approximation. This figure confirms that equi- 
librium theory predicts a fairly nontrivial temperature 
evolution of the anomaly location, which is itself a non- 
monotonic function of the temperature, with a crossover 
near T s=s 0.04. For T > 0.04 the anomaly shifts to larger 
density for increasing temperature, while the opposite 
behaviour is found for T < 0.04. 

The low temperature behaviour, T < 0.04, is read- 
ily understood on the basis of the equilibrium argument 
laid out above. When T decreases, the entropic penalty 



to find efficient packing of the particles contributes less 
to the free energy, and it is preferable to avoid the ener- 
getically costly particle overlaps up to larger volume frac- 
tions. Thus, from equilibrium considerations one would 
naturally predict qualitatively the low temperature be- 
haviour reported in Fig. [3J It reflects the physical intu- 
ition that thermal annealing should help the system to 
'solve the packing problem' more efficiently. 

Temperature has a second, competing effect which ex- 
plains the high temperature behaviour in Fig. [3j When 
T gets larger, particles acquire more kinetic energy, and 
they can more easily overlap. Thus the 'effective' particle 
diameter gets reduced by a simple thermal effect when T 
increases [H, [(| , and this effect must be compensated by 
an increase of the volume fraction. Therefore, tempera- 
ture has two opposite effects which combine to yield the 
behaviour reported in Fig. [3J and explain the crossover 
near T « 0.04. 

What happens if temperature is decreased further? We 
mentioned in the introduction that our equilibrium cal- 
culations eventually become inapplicable since a glass 
transition should take place when the structural relax- 
ation time gets very large. In a simulation or an experi- 
ment, there will be a temperature below which the system 
will eventually be frozen into some amorphous configura- 
tion . It is possible to describe this glass physics using 
liquid state theory, but this requires considerably more 
advanced techniques 
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Still, a number of conclusions naturally follow from the 
above equilibrium results. First, at low but finite temper- 
atures, the system is deep into the glass region. Although 
it cannot relax to equilibrium on accessible timescales, it 
attempts to do so, albeit at a very slow pace. Thus, the 
physical explanation given above for the very existence of 
the structural anomaly should still apply in this regime, 
which is the one explored in Refs. [l|, |4j. The structural 
anomaly should thus be 'frozen in' by the glass transi- 
tion, which involves no important change of the struc- 
ture. Second, we have shown above that the location of 
the maximum is ruled by two competing effects. In the 
glass phase, only the second one remains effective, as the 
system is not able to explore different configurations effi- 
ciently and find better particle packings. Thus, the evo- 
lution of the anomaly at low temperature should change 
again as the glass transition is crossed, and revert to the 
high temperature behaviour of Fig. [3j Overall, thus, we 
predict that the locus of the structural anomaly displays 
three regimes due to the competing effects of (i) thermal 
fluctuations reducing the effective particle diameter, (ii) 
the entropy / energy competition favouring more efficient 
packing at low temperatures, (iii) the intervening glass 
transition. In a (T, ip) phase diagram, the data should 
thus follow a nontrivial '5-shape'. This theoretical pre- 
diction is our main new result. 

We confront these theoretical considerations to numer- 
ical simulations in Fig. Q] We report the HNC prediction 
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FIG. 4: Location of the maximum of the first peak of g(r) in a 
temperature volume fraction phase diagram from HNC (filled 
circles), equilibrium simulations (open circles), nonequilib- 
rium glass configurations (open triangles). The glass tran- 
sition are defined as iso-relaxation timescales lines, with two 
different large values, from Ref. ||. The location of the max- 
imum follows a nontrivial <S-shape, as predicted theoretically. 



for the location of the structural anomaly from Fig. [3] 
and compare it to the outcome of molecular dynamics 
simulations at thermal equilibrium. To obtain these data 
points, we made use of the extensive set of data of Ref. Q , 
where a 50:50 binary mixture of harmonic spheres with 
size ratio 1.4 was studied using molecular dynamics. A 
mixture is needed in the simulation to avoid crystalliza- 
tion. When necessary, we performed additional simula- 
tions to explore a broader range of state points, using the 
same techniques. We use a large system size, N = 8000, 
where Unite size effects are negligliblc both at equilib- 
rium and at very low temperatures @, E3] ■ 

We find a very good qualitative agreement between the 
equilibrium behaviour predicted using HNC and the out- 
come of computer simulations at sufficiently high tem- 
peratures, T > 0.003, confirming that HNC performs 
quite well for harmonic spheres in this range of densities 
and temperatures. The nonmonotonic temperature be- 
haviour and the location of the crossover near T as 0.04 
are well reproduced. Quantitative agreement cannot be 
expected, as size polydispersity was not taken into ac- 
count in our theoretical calculations. 

We then report glass transition lines, defined from the 
temperature where the structural relaxation time, r, be- 
comes larger than some prescribed value. We report two 
lines corresponding to two large values, r = 10 3 and 
t = 10 4 , in the reduced time units used in molecular 
dynamics simulations Q. It is very difficult to maintain 
thermal equilibrium below these lines in computer sim- 
ulations. These glass lines meet the locus of structural 
anomaly near ip ss 0.8 and T as 0.0025, which shows that 
the maximum of g(r) cannot be observed in equilibrium 
conditions in computer simulations below T = 0.0025. 
This also suggests that both phenomena are not directly 



connected. 

At even lower temperatures, it becomes ambiguous to 
locate a unique maximum of g(r), as this should in prin- 
ciple depend on the way the glass configurations have 
been prepared. To illustrate this difficulty, and there- 
fore the noncquilibrium nature of the low temperature 
behaviour, we numerically prepared two glasses of har- 
monic spheres: 'Glass 1' is very rapidly quenched from a 
high temperature to well below the glass transition, while 
'Glass 2' is very slowly annealed across T g . For both 
glasses, we track the location of the g(r) maxima down 
to very low temperatures, see Fig.|U As discussed above, 
these data confirm the behaviour predicted at low tem- 
perature when no structural relaxation can occur. The 
location of the maximum shifts to lower volume fraction 
when temperature decreases. This is precisely the regime 
of temperature described in previous work [l|, [4| , and our 
results for 'Glass 1' agree with those earlier reports. In 
Fig. [4l we report the low temperature scaling behaviour 
deduced from considering the effect of thermal fluctua- 
tions alone, T ~ (ip — ip c ) 2 , where the exponent '2' is 
the same as in the pair potential |T]) and ip c denotes the 
location of the jamming transitions in the T = limit. 
These lines separate in the glass phase a region where the 
glass is dominated by repulsion at low volume fraction, 
from a region where particles are compressed at large 
volume fraction, and the system resembles a compressed 
emulsion. 

The noncquilibrium nature of the jamming transition 
is illustrated from the difference between our two glasses 
in Fig. [4j They both display a line of structural anomaly, 
but these two lines are distinct and define jamming tran- 
sitions at distinct densities, ip c = 0.65 and tp c = 0.662. 
The dependence of p c on the preparation protocol in the 
T = limit has been discussed before [ill, HI, H, H3 l- 



The data in Fig. [4] show that these considerations ap- 
ply to finite temperatures as well, and to the struc- 
tural anomaly in particular. Therefore, this figure illus- 
trates the fact that nonequilibrium jamming transition 
lines are as ambiguously defined as glass transition lines, 
since they both strongly depend upon operational defini- 
tions [12, H,[23. These two concepts can only become 
protocol-independent in theoretical calculations [H, [25| ■ 

The crossover between the very low temperatures in- 
vestigated in previous work, and the equilibrium regime 
studied here is also interesting. In the regime T as 
0.001 — 0.003, the system is an aging glass at large den- 
sity, but a fluid at low density. For the fluid state points 
the first peak of g(r) increases with ip and the extrapo- 
lated maximum is at very large tp. For the aging glasses, 
however, the peak decreases with ip as the glass maximum 
should follow the low temperature trend. Thus, these two 
tendencies produce a maximum when the glass transition 
is crossed, and the glass and jamming lines in Fig. [U 'kiss' 
in this temperature regime. Of course, the exact location 
of the maximum depends on the aging time. This effect 
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was not reported before [H, Q , as the structure docs not 
age very much at the very low temperatures studied in 
these earlier reports, T/T g « 1/1000 - 1/10. 

It is interesting to remark the close similarity between 
the effect of glass annealing in nonequilibrium situations, 
and the effect of decreasing the temperature at equilib- 
rium: they both shift the maximum of g(r) to larger <p for 
the same physical reason. This correspondence empha- 
sizes once more that the physics of the structural anomaly 
described in this paper is in fact deeply rooted in equilib- 
rium concepts. The T — > limit is no exception, but it 
makes the distinction between repulsive and compressed 
regimes sharper, as the energy of the ground state of the 
system is either or positive Q. 

To conclude, we have used a simple equilibrium theo- 
retical approach to study the existence and location of a 
structural anomaly reported in recent experiments in soft 
repulsive spheres in nonequilibrium conditions. We have 
suggested that the structural anomaly reported in these 
articles can be explained using purely equilibrium con- 
cepts, and several aspects of its behaviour are indeed well 
described using liquid state theory. The very existence 
of an anomaly relies on a standard competition between 
entropy and energy in the free energy of these very soft 
particles, and we have suggested that the anomaloy be- 
longs to a much broader class of similar anomalies studied 
in the literature. In the present case, it directly results 
from the extreme softness of particles described by pair 
potentials such as Eq. (JTJ. 

For standard models describing atomic liquids, such 
as Lennard- Jones potentials, the repulsive core is much 
steeper, and they are thus governed by a different physics 
upon compression. Indeed, the large density limit of both 
types of systems is very different, as the harmonic spheres 
behave as an ideal gas in this limit [l?! • In the context of 
glassy physics, a recent quantitative comparison between 
harmonic spheres and Lennard- Jones models of the glass 
transition has also confirmed the different nature of their 
glassy dynamics and of its evolution with density [28| . 

Additionally, we have discussed three competing ef- 
fects of the temperature which control the location of the 
anomaly in the phase diagram, and predicted a nontrivial 
'iS-shapc', confirmed by numerical simulations. We have 
finally shown that the anomalous structural behaviour of 
soft particles discussed in Refs. [l], 0, d, Q m fact occurs 
over a broad, continuous range of densities and temper- 
atures, as it is directly affected by sample preparation 
effects. 

In future work [2BI ] , we plan to extend the present liquid 
state calculations to describe analytically the properties 
of dense systems of soft particles down to very low tem- 
peratures, where ergodicity breaking must explicitely be 
taken into account. 
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